Celem przedmiotowego ćwiczenia jest zapoznanie się studentów z pracą z danymi pochodzącymi z bazy GIOś oraz z danymi meteorologicznymi. Aby prawidłowo wykonac zlecone nam zadanie musielismy zapoznac sie z formatem danych, poznać narzędzia które ułatwią nam pracę oraz ich przetwarzanie. Musieliśmy dokonać prognozy wybranego przez nasz zanieczyszczenia powietrza na dany rok, oraz dokonać analizy w celu określenia największej korealcji danych aby z łatwością móc odczytać panujące trendy w wybranej stacji obserwacyjnej. Ostatnim aspektem naszej pracy, było napisanie sprawozdania technicznego z wykonanengo ćwiczenia przy pomocy RMarkdown, w celu udokumentowania naszej pracy w języku R.
R - interpretowany jezyk programowania oraz srodowisko do obliczen statystycznych. Stosowany jest w analizie szeroko rozumianych danych środowiskowych i przestrzennych oraz ich wizualizacji. Podobny jest do języka i środowiska S stworzonego w Bell Laboratories przez Johna Chambersa i jego współpracowników. R jako implementacja języka S została stworzona przez Roberta Gentlemana i Rossa Ihakę na uniwersytecie w Auckland. Nadaje się on świetnie do interaktywnej pracy z danymi, ponieważ połączono w nim wybrane cechy języków funkcyjnych oraz obiektowych.
R Markdown - jest formatem pliku stworzonym do sporzadzania dynamicznych dokumentow z wykorzystaniem R. Plik typu Markdown jest pisany w specyficzny dla siebie sposob, który zaklada bardzo latwa edycje tekstu oraz implementowanie w nim fragmentow kodu (chunki zawierajace kod z poleceniami w jezyku R). R Markdown jest bardzo wygodna metoda formatowania plikow HTML, PDF i dokumentow MS Word.
Baza GIOŚ - (Główny Inspektorat Ochrony Środowiska) - to zbiór publicznie dostępnych danych zebranych przez organ jakim jest GIOŚ.
Jednym z najistotniejszych zadań realizowanych przez Inspekcję Ochrony Środowiska jest prowadzenie badań i ocen stanu środowiska, w tym monitoringu jakości powietrza. Zadanie to jest wykonywane w ramach Państwowego Monitoringu Środowiska (PMŚ), którego program jest opracowywany przez Głównego Inspektora Ochrony Środowiska i zatwierdzany przez Ministra Środowiska. W oparciu o krajowy program PMŚ opracowywane są wojewódzkie programy PMŚ zatwierdzane przez Głównego Inspektora Ochrony Środowiska.
Monitoring jakości powietrza obejmuje zadania związane z badaniem i oceną stanu zanieczyszczenia powietrza, w tym pomiary i oceny jakości powietrza w strefach, monitoring tła miejskiego pod kątem WWA, pomiary stanu zanieczyszczenia powietrza pyłem PM2,5 dla potrzeb monitorowania procesu osiągania krajowego celu redukcji narażenia, pomiary stanu zanieczyszczenia powietrza metalami ciężkimi i WWA oraz rtęcią w stanie gazowym na stacjach monitoringu tła regionalnego, pomiary składu chemicznego pyłu PM2,5, monitoring prekursorów ozonu; programy badawcze dotyczące zjawisk globalnych i kontynentalnych wynikające z podpisanych przez Polskę konwencji ekologicznych.
Ok. 90% pomiarów jakości powietrza wykonywanych w ramach PMŚ oraz roczne i pięcioletnie oceny jakości powietrza w strefach są wykonywane przez GIOŚ i regionalne wydziały. Na zlecenie Głównego Inspektoratu Ochrony Środowiska są realizowane krajowe programy monitoringu jakości powietrza, GIOŚ jednocześnie nadzoruje i koordynuje wykonywanie programu badań i ocen jakości powietrza określonego w krajowym i wojewódzkich programach Państwowego Monitoringu Środowiska.
Las Losowy (ang. Random Forest) - jedna z metod ensamblingu uczenia maszynowego. Jest to rozszerzenie metody drzew decyzyjnych, która polega na budowaniu modelu w oparciu o zadawanie i dzielenie zbioru według warunków nakładanych na dane zmienne. I chociaż lasy losowe nadają się lepiej do prognozowania danych dyskretnych bądź nominalnych (metoda klasyfikacyjna), to również w naszym przypadku powinna przynieść satysfakcjonujące efekty.
[Źródlo] Materiał Zajęciowy
Aby rozpocząc pracę z naszymi danymi, musieliśmy je pozyskać, w tym celu skorzystaliśy z środowiska dołączonego przez Pana Doktora, jako materiał zajęciowy “projekt_2.RData”. Aby zaimportować dane skorzystalićmy z poniższego polecenia:
load(file = "projekt_2.RData")Naszą pracę z danymi postanowiliśmy rozpocząć od usnięcia polskich znaków w kolumnach z których w dalszczej części projektu będziemy korzystać:
colnames(gios_inv)[c(7,13)] <- c("Data.zamkniecia", "Miejscowosc")
colnames(gios_inv)[15:16] <- c("Szerokosc", "Dlugosc")W celu wykonania naszego ćwiczenai postanowiliśmy skorzystać z stacji znajdującej się północnej cześci Częstochowy, w dzielnicy “Tysiąclecia” na ulicy. Krzysztofa Kamila Baczyńskiego. W tym celu skorzystaliśmy z stacji SlCzestoBacz.
gios_inv %>% filter(Kod.stacji == "SlCzestoBacz") -> czesto_lok
czesto_lok_info <- paste(paste(
czesto_lok$Kod.stacji,
paste("Miejscowość:", czesto_lok$Miejscowosc),
paste("Data uruchomienia:", czesto_lok$Data.uruchomienia),
paste("Data zamknięcia:", czesto_lok$Data.zamkniecia),
paste("Typ stacji:", czesto_lok$Typ.stacji),
paste("Typ obszaru:", czesto_lok$Typ.obszaru),
paste("Współrzędne:", czesto_lok$Dlugosc,"E, ", czesto_lok$Szerokosc,"N"),
sep = "<br/>"
))
leaflet() %>%
addTiles() %>%
addMarkers(data=czesto_lok,
lng= ~ Dlugosc,
lat= ~ Szerokosc,
popup = czesto_lok_info)Lokalizację geograficzną naszej stacji jakości powietrza sprawdziliśmy nieco wcześniej aby być pewnymi że w jej poblizu znajduje się stacja meteorologiczna z której będziemy mogli skorzystać dla prawidłowej analizy. Jak widzimy na poniższej mapie, odległość dzieląca nasze stacje nie przekracza 4km, dzięki czemu dane będą bardzo obiektywne.
noaa_isd <- getMeta(end.year="current", lon=19.130111, lat=50.836389, returnMap=T)
noaa_isdPo znalezieniu interesujacej nas stacji, postnaowiliśy zaimportować interesujące nas dane do zmiennej, tak aby w dalszej częsci naszej pracy nad projektem, z łatwością móc z nich korzystać, bez wracanie do tego momentu przygotowywania danych.
czestochowa_met <- importNOAA(code="125500-99999", year=2000:2020)Postanowiliśmy dokonać analizy ciagłości obserwacji, dla wybranej przez nas stacji, aby tego dokonać skorzystaliśmy z poniższego fragmentu kodu:
# Dodanie informacji o roku do ramki danych
czestochowa_met$years <- format(czestochowa_met$date, "%Y")
# Policzenie ilości pomiarów dla danego roku dla stacji Okęcie
czestochowa_met %>%
group_by(years) %>%
summarise(liczba_pomiarow = n() - sum(is.na(wd))) -> czestochowa_met_summary
czestochowa_met_summary <- czestochowa_met_summary %>%
as_tibble() %>%
mutate(years = as.integer(years))
czy_pomiary_czestochowa_met <- data.frame(ymin = -Inf,
ymax = Inf,
xmin = c(1999 , 2021),
xmax = c(1999.5 , 2021.5))
# Utworzenie wykresu
ggplot(czestochowa_met_summary, aes(x=years, y=liczba_pomiarow, color=liczba_pomiarow)) +
geom_point(size=3, show.legend = FALSE) +
scale_y_continuous(name = 'Liczba rekordów pomiarowych',
limits = c(0, 9000),
expand = c(0, 0)) +
scale_x_continuous(name = 'Rok',
limits = c(1998, 2022),
expand = c(0, 0),
breaks = seq(2000, 2020, by = 1)) +
geom_rect(data = czy_pomiary_czestochowa_met,
aes(ymin = ymin, ymax = ymax,
xmin = xmin, xmax = xmax), fill = 'red',
alpha = 0.15, inherit.aes = FALSE) +
theme_minimal() +
theme(plot.margin = unit(c(.5,.5,.5,.5), "cm")) +
labs(title = 'Liczba pomiarów dla stacji pomiarowej CZESTOCHOWA (CZESTOCHOWA) w ISD NOAA',
subtitle = paste0('Stan na ',format(Sys.time(),'%d %B, %Y'))) Stacja Częstochowa działa od 1940 roku, w badanym przez nas okresie czasu - to jest od 2000 - 2020 roku, stacja nie posiadała przerwy w dostarczaniu danych meteorologicznych. dodatkowo z powyższego wykresu możemy zauważyć, że najlepszą dokładnością danych charakteryzuje się okres od 2017-2019, gdyż to właśnie w tym okresie wykonywanych było najwięcej pomiarów, oraz zebrane dane miały charakter ciągły.
Analizę naszych danych rozpoczeliśmy od wyselkecjonowania interesujących nas danych dotyczących zanieczyszczenia pyłami PM10 z wybranej przez nas stacji. Pozwoli nam to znacząco ograniczyć czas wykonywanych poleceń, oraz pozwoli skupić się na dokładniejszej analizie wybranego zagadnienia.
czestochowa_PM10 <- PM10_1h %>%
filter(kod=="SlCzestCzes_baczy" | kod=="SlCzestoBacz") %>%
select(date, obs) Aby nasze dane prawidłowo z sobą współgrały musieliśy wykonać szereg poleceń mających na celu dodanie godziny z uwzględnieniem stefy czasowej, dokonać selekcji wybranych danych a następnie połaczyć je w jeden plik z danymi. W tym celu skrozystaliśmy z poleceń :
## Musimy dodać godzinę w danych meteo, aby uwzględnić różne strefy czasowe danych
czestochowa_met$date <- czestochowa_met$date + 3600
## selekcja wybranych zmiennych
czestochowa_met <- czestochowa_met[c(3, 7:14)]
## Łączenie danych meteorologicznych z danymi jakości powietrza
dane <- inner_join(czestochowa_met, czestochowa_PM10, by = "date")Nasze dane dla lepszego zobrazowania postanowiliśmy przedstawić w formie tabeli danych aby każdy mógł odczytać interesujące go dane
datatable(dane)W rozpatrywanym okresie występuje wysoka kompletnosć danych. W dalszej czesci naszej analizy posluzymy sie wynikami w formie graficznej, gdyz przy takiej ilosci danych, taka forma bedzie bardziej czytelna.
Wykorzystując to metodę niezbędne jest posiadanie kompletnych wierszy obserwacji, a także tylko tych zmiennych, które wniosą informację dodaną do naszego modelu. Aby upewnić się w którym roku posiadanmy najlepszą kompletność danych, skorzystaliśmy z poniższego polcenia, jednocześnie wybierając te które charakteryzowały się najlepszą kompletnością.
# Sprawdzenie braków danych i poszczególnych kolumn
summary(dane)## date ws wd air_temp
## Min. :2006-01-01 00:00:00 Min. : 0.00 Min. : 0.0 Min. :-25.70
## 1st Qu.:2009-04-01 20:00:00 1st Qu.: 2.00 1st Qu.:140.0 1st Qu.: 2.20
## Median :2012-07-01 17:00:00 Median : 2.00 Median :210.0 Median : 9.00
## Mean :2012-07-01 20:36:55 Mean : 2.48 Mean :201.6 Mean : 9.06
## 3rd Qu.:2015-10-02 00:00:00 3rd Qu.: 3.00 3rd Qu.:280.0 3rd Qu.: 16.00
## Max. :2018-12-31 23:00:00 Max. :13.00 Max. :360.0 Max. : 36.50
## NA's :42861 NA's :43668 NA's :42216
## atmos_pres visibility dew_point RH
## Min. : 976 Min. : 0 Min. :-28.60 Min. : 13.03
## 1st Qu.:1012 1st Qu.: 5000 1st Qu.: -0.50 1st Qu.: 65.10
## Median :1017 Median :10000 Median : 5.00 Median : 81.98
## Mean :1017 Mean :12995 Mean : 4.61 Mean : 77.00
## 3rd Qu.:1022 3rd Qu.:20000 3rd Qu.: 10.50 3rd Qu.: 92.38
## Max. :1051 Max. :50000 Max. : 21.90 Max. :100.00
## NA's :42250 NA's :42316 NA's :42250 NA's :42250
## ceil_hgt obs
## Min. : 15 Min. : 0.00
## 1st Qu.: 750 1st Qu.: 15.36
## Median : 1200 Median : 25.43
## Mean : 8841 Mean : 35.81
## 3rd Qu.:22000 3rd Qu.: 43.00
## Max. :22000 Max. :705.91
## NA's :92012 NA's :13122
# Selekcja tylko kompletnych wierszy
dane_rf <- dane_rf[complete.cases(dane_rf),]W celu przeprowadzenia dokładniejszej analizy postanowiliśmy przeprowadzić analizę dla więcej niż jednego modelu. W tym celu postanowiliśmy utworzyć nowe zmienne naszych danych różniące się parametrami, aby otrzymać różny punkt odniesienia, na którym będziemy mogli bazować. Zdecydowaliśmy się na utworzenie 5 zmiennych gdyż taka ilość nie powodowała błędów wykonywanego kodu, oraz czas pracy wykonywanej przez system operacyjny był stosunkowo szybki.
dane_rf %>% mutate(years = year(date)) -> dane_rf
#dodanie kolumn
dane_rf$unix_date <- as.numeric(as.POSIXct(dane_rf$date))
dane_rf %>% mutate(years = year(date)) -> dane_rf
dane1 <- dane_rf %>%
mutate(week = week(date))
dane2 <- dane1 %>%
mutate(hour = hour(date))
dane3 <- dane2 %>%
mutate(weekday = wday(date))
dane4 <-dane3 %>%
mutate(month = month(date))Wykreślmy wykres, który poznaliśmy w czasie wykonywania Projektu numer jeden, który ma na celu zobrazowanie ilości rekordów w danym roku:
# Policzenie ilości pomiarów w poszczególnych latach
dane_rf$years <- format(dane_rf$date, "%Y")
dane_rf %>%
group_by(years) %>%
summarise(liczba_pomiarow = n()) -> dane_rf_summary
ggplot(dane_rf_summary, aes(x=years, y=liczba_pomiarow, color=liczba_pomiarow)) +
geom_point(size=3, show.legend = FALSE) +
theme_minimal() +
coord_flip()Zgodnie z przykładowym projektem II, dołączonym jako materiał zajęciowy, dzielimy wszystkie zbiory na dwa podzbiory. Jeden z nich, tzw. “zbiór treningowy”, będzie obejmować dane, na których nauczymy model przewidywać stężenia PM10 na podstawie danych meteorologicznych. Nie wybierajmy zbyt wiele rekordów, bo zwiększy to wykładniczo czas obliczeniowy. Zdecydujmy się na 2 - 3lata najbardziej kompletnych danych. My do zbioru treningowego postanowiliśmy skorzystać z lat 2016 - 2018. Do zbioru, testowego, czyli takiego, na którym będziemy weryfikować nasz model, wybieramy rok 2018. Właśnie dla tego przedziału skorzystaliśmy dla każdej z naszych zmiennych.
# Filtrujemy dane do zbioru treningowego i testowego
dane_rf %>%
filter(years>2015 & years<2018) %>%
select(-date, -years) -> dane_rf_train
dane_rf %>%
filter(years==2018) %>%
select(-years) -> dane_rf_test
dane1 %>%
filter(years>2015 & years<2018) %>%
select(-date, -years) -> dane1_rf_train
dane1 %>%
filter(years==2018) %>%
select(-years) -> dane1_rf_test
dane2 %>%
filter(years>2015 & years<2018) %>%
select(-date, -years) -> dane2_rf_train
dane2 %>%
filter(years==2018) %>%
select(-years) -> dane2_rf_test
dane3 %>%
filter(years>2015 & years<2018) %>%
select(-date, -years) -> dane3_rf_train
dane3 %>%
filter(years==2018) %>%
select(-years) -> dane3_rf_test
dane4 %>%
filter(years>2015 & years<2018) %>%
select(-date, -years) -> dane4_rf_train
dane4 %>%
filter(years==2018) %>%
select(-years) -> dane4_rf_testAby w przyszłosći nie pojawił się problem spójności danych, postanowiliśmy zapisać wyselekcjonowane dane jako nowa zmienna. W dalszej analizie bedziemy bazować na, kopii tak aby nie uszkodzić pliku oryginalnego.
# Zapisanie kopii obiektu, tak aby nie zmieniać nazw kolumn w oryginalnym zestawie danych
dane_rfc <- dane_rf_train
dane1_rfc <- dane1_rf_train
dane2_rfc <- dane2_rf_train
dane3_rfc <- dane3_rf_train
dane4_rfc <- dane4_rf_trainAby sprawdzić korelacje liniowe między naszymi zmiennymi posłużymy się korelacją Pearsona - tą samą z której krozystaliśmy w Projekcie I.
Macierz korelacji liniowej Pearsona (macierz współczynników określających poziom zależności liniowej między zmiennymi losowymi), pozwala wstępnie sprawdzić występowanie zależności między analizowanymi zmiennymi. Mówiąc prosto, współczynnik ten może przyjmować wartości od -1 do 1. Im jest mniejszy, tym silniejsza zależność ujemna między parametrami, a im jest bliższy 1, tym silniejsza zależność dodatnia. Dodatkowo możemy sprawdzić, czy otrzymane wartości są istotne statystycznie (przyjmiemy poziom istotności równy 0.05).
Aby mówić w ogóle o wykorzystaniu modelu, musimy mieć solidne podstawy, aby sądzić, że związek między zmiennymi jest silny i jedną z nich (PM10 - zmienna objaśniana) da się przedstawić jako funkcja pozostałych zmiennych (zmienne objaśniające).
# Definiowanie indywidualnej palety kolorów w zależności od współczynnika korelacji (kodowanie hex)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
# Policzenie odpowiednich współczynników za pomocą funkcji 'corr'
corr <- rcorr(as.matrix(dane_rfc))
# Zdefiniowanie wektora ze współczynnikami korelacji
corr_r <- corr$r
# zdefiniowanie wektora z wartościami p-value
corr_p <- corr$P
# Wykreślenie macierzy korelacji
corrplot(corr_r, method = "color", col = col(200),
type = "upper", order = "hclust", addCoef.col = "black",
diag = FALSE,
tl.col = "black",
tl.srt = 45,
p.mat = corr_p,
sig.level = 0.05)Jak widzimy, korelacje liniowe Pearsona między stężeniem pyłu PM10 a innymi zmiennymi są jedynie na umiarkowanym poziomie, ale wszystkie są statystycznie istotne. Ta wiedza jest dla nas wystarczająca, aby być pewnym że na naszych danych możemy zbudować model, jednocześnie nie musimy przejmować się, że jakaś korealacja okaże się nieistotna statystycznie - model prawdopodobnie sam ją odrzuci na etapie uczenia.
Do uruchomienia lasu losowego wykorzystano pakiet caret, który jest bardzo bogatym pakietem do modelowania statystycznego w R. Aby oszczędzić nieco czasu zgodnie z propozycją zawartą w konspekcie uruchomimy także obliczenie równoległe, które wykona się wykorzystując dostępne rdzenie procesora. Skorzystaliśmy z takiej metody gdyż utworzenie 5 modeli, było bardzo czasochłonne na naszym sprzęcie, a ta metoda pozwoliła odrobinę oszczędzić nam czasu, oraz ograniczyłą zużycie zasobów naszego sprzętu.
# Ustawiamy metodę tzw. "kroswalidacji", czyli samo-testowania się modelu w trakcie jego budowy.
cross.walid <- trainControl(method = "cv",
number = 5,
allowParallel = T,
returnResamp = 'all')
# Zadajemy, aby w modelu były obecne wszystkie zmienne oprócz objaśnianej (PM10)
tunegrid <- expand.grid(.mtry = c(1:10))
tunegrid1 <- expand.grid(.mtry = c(1:11))
tunegrid2 <- expand.grid(.mtry = c(1:12))
tunegrid3 <- expand.grid(.mtry = c(1:13))
tunegrid4 <- expand.grid(.mtry = c(1:14))
tunegrid5 <- expand.grid(.mtry = c(1:15))
# Uruchamiamy klaster obliczeniowy z wszystkich dostępnych rdzeni - 1 (zwyczajowo zostawia się jeden, aby system swobodnie działał)
cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)
# Budujemy (trenujemy) model lasu losowego
model_rf <- train(obs ~ .,
data = dane_rf_train,
method = "rf",
replace = T,
ntree = 200,
metric = 'RMSE',
tuneGrid = tunegrid,
trControl = cross.walid)
model1_rf <- train(obs ~ .,
data = dane1_rf_train,
method = "rf",
replace = T,
ntree = 200,
metric = 'RMSE',
tuneGrid = tunegrid1,
trControl = cross.walid)
model2_rf <- train(obs ~ .,
data = dane2_rf_train,
method = "rf",
replace = T,
ntree = 200,
metric = 'RMSE',
tuneGrid = tunegrid2,
trControl = cross.walid)
model3_rf <- train(obs ~ .,
data = dane3_rf_train,
method = "rf",
replace = T,
ntree = 200,
metric = 'RMSE',
tuneGrid = tunegrid3,
trControl = cross.walid)
model4_rf <- train(obs ~ .,
data = dane4_rf_train,
method = "rf",
replace = T,
ntree = 200,
metric = 'RMSE',
tuneGrid = tunegrid4,
trControl = cross.walid)
# Zatrzymujemy klaster
stopCluster(cluster)
registerDoSEQ() Z poniższych wykresów możemy z łatwością zobaczyć zmiany w testowaniu się modelu.
#Model 9 zmiennych
plot(model_rf)#Model 10 zmiennych
plot(model1_rf)#Model 11 zmiennych
plot(model2_rf)#Model 12 zmiennych
plot(model3_rf)#Model 13 zmiennych
plot(model4_rf)Jednak aby dokładnie odczytać wartość mtry, RMSE, Rsquared oraz MAE, skorzystamy z poniższych wyników .
#Model 9 zmiennych
model_rf## Random Forest
##
## 11184 samples
## 8 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8947, 8947, 8947, 8947, 8948
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 1 22.14724 0.8040095 12.94263
## 2 20.82950 0.8196639 12.17065
## 3 20.58272 0.8228228 12.03239
## 4 20.50626 0.8232116 11.96935
## 5 20.49883 0.8227717 11.95867
## 6 20.45653 0.8231178 11.94832
## 7 20.59064 0.8205505 12.02360
## 8 20.75866 0.8173084 12.09947
## 9 20.70904 0.8182752 12.08131
## 10 20.70458 0.8183340 12.07470
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 6.
#Model 10 zmiennych
model1_rf## Random Forest
##
## 11184 samples
## 10 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8948, 8946, 8948, 8947, 8947
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 1 19.53636 0.8493156 11.33375
## 2 18.16269 0.8648520 10.50227
## 3 18.06996 0.8647392 10.41039
## 4 18.13147 0.8631444 10.43246
## 5 18.17947 0.8616916 10.44163
## 6 18.22309 0.8609007 10.45905
## 7 18.23165 0.8605247 10.51912
## 8 18.37433 0.8582054 10.55277
## 9 18.35812 0.8581134 10.60079
## 10 18.58768 0.8544083 10.69067
## 11 18.58530 0.8542627 10.69724
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 3.
#Model 11 zmiennych
model2_rf## Random Forest
##
## 11184 samples
## 11 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8947, 8948, 8948, 8946, 8947
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 1 19.64163 0.8515474 11.46933
## 2 17.82246 0.8707086 10.30115
## 3 17.60785 0.8722647 10.20005
## 4 17.59252 0.8721022 10.19026
## 5 17.64774 0.8705475 10.22497
## 6 17.70192 0.8695952 10.21369
## 7 17.66486 0.8699411 10.23112
## 8 17.84399 0.8669006 10.30405
## 9 17.77658 0.8679943 10.31162
## 10 17.86031 0.8660553 10.37553
## 11 18.05088 0.8630838 10.44521
## 12 17.96967 0.8642814 10.41103
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
#Model 12 zmiennych
model3_rf## Random Forest
##
## 11184 samples
## 12 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8948, 8947, 8946, 8947, 8948
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 1 19.32859 0.8604909 11.25531
## 2 17.52587 0.8771690 10.18301
## 3 17.30842 0.8778948 10.05658
## 4 17.19356 0.8783876 10.00094
## 5 17.31842 0.8760826 10.04116
## 6 17.32375 0.8754202 10.06966
## 7 17.32627 0.8748956 10.09337
## 8 17.50574 0.8720442 10.14486
## 9 17.53227 0.8714417 10.15545
## 10 17.52422 0.8713014 10.19491
## 11 17.62599 0.8694167 10.24135
## 12 17.69618 0.8681833 10.28774
## 13 17.65027 0.8685643 10.28715
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
#Model 13 zmiennych
model4_rf## Random Forest
##
## 11184 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 8947, 8948, 8948, 8947, 8946
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 1 19.81016 0.8546864 11.395945
## 2 17.61938 0.8761540 10.080288
## 3 17.44280 0.8760401 9.939129
## 4 17.35046 0.8767173 9.907598
## 5 17.40871 0.8753637 9.948455
## 6 17.52787 0.8730406 10.006474
## 7 17.51726 0.8731641 10.017711
## 8 17.58470 0.8716453 10.070079
## 9 17.66427 0.8703834 10.105871
## 10 17.65625 0.8705696 10.129317
## 11 17.79102 0.8683707 10.196462
## 12 17.79616 0.8682362 10.220303
## 13 17.95914 0.8654848 10.306924
## 14 17.94380 0.8657578 10.313816
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
Wykorzystując przed chwilą storzone modele możemy sprawdzić jak zachowują się one na danych, które nie były wykorzystane do jego uczenia. Pobierane zostaną dane meteorologiczne dla roku 2018 dla tej samej stacji w Częstochowie, a następnie wykorzystany zostanie utworzony model do oszacowania stężenia pyłu PM10 w tym roku.
dane_rf_test$mod <- predict(model_rf, dane_rf_test)
dane1_rf_test$mod <- predict(model1_rf, dane1_rf_test)
dane2_rf_test$mod <- predict(model2_rf, dane2_rf_test)
dane3_rf_test$mod <- predict(model3_rf, dane3_rf_test)
dane4_rf_test$mod <- predict(model4_rf, dane4_rf_test)Żeby sprawdzić, jak dobrze nasz model odwzorowuje rzeczywistość musimy mieć jakąś wartość odniesienia. Sporządźmy dwa wykresy obrazujące nam zależność wartości stężenia PM10 zaobserwowanego na stacji jakości powietrza od zamodelowanej jego wartości. Pierwszy z nich to wykres rozrzutu.
# Wykres rozrzutu
# 9 zmiennych
ggplot(data=dane_rf_test, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
theme_minimal()# 10 zmiennych
ggplot(data=dane1_rf_test, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
theme_minimal()# 11 zmiennych
ggplot(data=dane2_rf_test, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
theme_minimal()# 12 zmiennych
ggplot(data=dane3_rf_test, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
theme_minimal()# 13 zmiennych
ggplot(data=dane4_rf_test, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
theme_minimal()Wykres rozrzutu pokazuje jak mocno wyniki prognoz odchylone są od pomiarów stężeń PM10. Szara przerywana linia obrazuje model idealny (w sytuacji, gdy prognoza równa jest pomiarowi, czyli x=y), a niebieska – linię trendu uzyskaną z naszych danych, o równaniu wyświetlonym na wykresie.
Według naszych danych, prognoza jest zaniżona w stosunku do pomiarów stężenia PM10 na stacji w Częstochowie, co jest szczególnie widoczne w przypadku najwyższych stężeń. Może to wynikać z faktu, że warunki meteorologiczne nie są jedynymi czynnikami wpływającymi na stężenie zanieczyszczenia na danym obszarze. Pominęliśmy całkowicie kwestię emisji, która, w przypadku sytuacji wyjątkowych i chwilowych, może drastycznie wpłynąć na mierzone stężenia zanieczyszczeń.
Uznaliśmy że porównanie tych wykresów kiedy znajdują się jeden pod drugim, jest bardzo nie praktyczne. W tym celu postanowiliśy przedstawić je bardziej czytelnie. W tym celu postanowiliśy opisać każdy z wykresóW zgodnie z jego modelem, a następnie zrobić jeden obraz zbiorczy, w tym celu skrozystalisy z poniższych poleceń:
bind_rows(dane_rf_test %>%
mutate(typ = 'Model A [+unix_date]'),
dane1_rf_test %>%
mutate(typ = "Model B [+week]"),
dane2_rf_test %>%
mutate(typ = "Model C [+hour]"),
dane3_rf_test %>%
mutate(typ = "Model D [+weekday]"),
dane4_rf_test %>%
mutate(typ = "Model E [+month]")) -> new_data
ggplot(data=new_data, aes(x=obs, y=mod))+
geom_point(alpha=0.6, color="purple") +
geom_smooth(method="lm", formula = y~x-1) +
geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0,
formula = y~x-1, parse = TRUE, size = 5) +
xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nd stacji w Czestochowie w roku 2018") +
theme_minimal() +
facet_wrap(~typ)Dzięki takiej formie prezentacji możemy, z łatwością porównać utworzone modele, odczytać interesujace nas parametry, oraz określić który model charakteryzuje się najlepszą dokładnoscią.
Wykresy dodatkowo dostarcza informacji o równaniu regresji to znaczy o równaniu linii dopasowanej do danych przedstawionych na wykresie. Równanie dla najlepszego wyniku naszego doświadczenia przedstawia się następująco: y = 0, 841x.
Dla danych powyższych wykresu wyliczono też wartość współczynnika determinacji oznaczonego jako R2. Współczynnik ten w przypadku naszych wykresów jest równy = 0,81 (w kilku przypadkach). Zgodnie z przedstawioną poniżej tabelą można wyciagnąć wnioski że rozrzut naszej prognozy jest dobrze dopasowany do pomiaru idealnego.
R2## # A tibble: 5 x 2
## X1 X2
## <chr> <chr>
## 1 0.0 – 0.5 dopasowanie niezadowalające
## 2 0.5 – 0.6 dopasowanie słabe
## 3 0.6 – 0.8 dopasowanie zadowalające
## 4 0.8 – 0.9 dopasowanie dobre
## 5 0.9 – 1.0 dopasowanie bardzo dobre
Ostantim aspektem naszej pracy jest utworzenie wykresu trendu prognozy dla obserwacji steżeń PM10. Postanowiliśmy utworzyć jeden wykres, który dzięki wysokiej dokładności, powinen być najbardziej dokłądną reprezentacją danych. Skorzystalśmy z zmiennej dane4_rf_test, gdyż to właśnie dla tych danych wyrkes rozrzuty charakteryzował się najlepszą dokładnością.
# Porównanie serii danych w danym fragmencie roku
ggplot(data = dane4_rf_test) +
geom_path(aes(x = date, y = obs, color = "Obserwacje"), size = 1.5) +
geom_path(aes(x = date, y = mod, color = "Model"), size = 1.5, alpha=0.6) +
scale_x_datetime(limits = ymd_h(c("2018-01-01 00", "2018-01-31 23"))) +
scale_y_continuous(limits=c(0,250)) +
scale_color_manual(values=c("Obserwacje"="grey","Model"="red"))+
xlab("Data") +
ylab("Stężenie PM10 [µg/m3]") +
ggtitle("Wykres trendu prognozy od obserwacji dla stężeń PM10 na stacji w Częstochowie w roku 2018") +
theme_minimal() +
theme(legend.position="top",legend.title=element_blank())Wykres trendu prognoz (linia czerwona) od obserwacji stężeń PM10 (linia szara) pozwala na ogólną ocenę zgodności tych dwóch serii danych. Możemy zauważyć, że wniosek wyciągnięty na podstawie poprzedniego wykresu jest potwierdzony – najwyższe wartości stężeń są zaniżane w naszym modelu. Ogólny trend jest zachowany, chociaż również dla najniższych stężeń zdarzają się odchylenia.
Oprócz grafiki wizualizacji oceny dokąłdności modelu zawsze warto przyjrzeć się paramtrom statystycznym oceny dokładności prognoz.
# 9 zmiennych
modStats(mydata = dane_rf_test,
mod = "mod",
obs = "obs",
type = "season") %>%
knitr::kable(digits = 2)| season | n | FAC2 | MB | MGE | NMB | NMGE | RMSE | r | COE | IOA | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | spring (MAM) | 1413 | 0.84 | -9.37 | 17.77 | -0.18 | 0.34 | 28.85 | 0.79 | 0.41 | 0.70 |
| 3 | summer (JJA) | 1185 | 0.85 | -3.62 | 8.51 | -0.15 | 0.35 | 12.07 | 0.42 | 0.10 | 0.55 |
| 1 | autumn (SON) | 1276 | 0.82 | -5.69 | 16.22 | -0.13 | 0.37 | 24.08 | 0.56 | 0.24 | 0.62 |
| 4 | winter (DJF) | 1945 | 0.85 | 0.18 | 20.02 | 0.00 | 0.37 | 31.52 | 0.70 | 0.39 | 0.69 |
# 10 zmiennych
modStats(mydata = dane1_rf_test,
mod = "mod",
obs = "obs",
type = "season") %>%
knitr::kable(digits = 2)| season | n | FAC2 | MB | MGE | NMB | NMGE | RMSE | r | COE | IOA | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | spring (MAM) | 1413 | 0.88 | -5.80 | 17.28 | -0.11 | 0.33 | 29.16 | 0.77 | 0.42 | 0.71 |
| 3 | summer (JJA) | 1185 | 0.81 | 4.22 | 9.73 | 0.17 | 0.40 | 12.70 | 0.37 | -0.03 | 0.49 |
| 1 | autumn (SON) | 1276 | 0.86 | -0.55 | 15.80 | -0.01 | 0.36 | 22.25 | 0.60 | 0.26 | 0.63 |
| 4 | winter (DJF) | 1945 | 0.79 | 6.60 | 22.02 | 0.12 | 0.40 | 31.88 | 0.70 | 0.33 | 0.66 |
# 11 zmiennych
modStats(mydata = dane2_rf_test,
mod = "mod",
obs = "obs",
type = "season") %>%
knitr::kable(digits = 2)| season | n | FAC2 | MB | MGE | NMB | NMGE | RMSE | r | COE | IOA | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | spring (MAM) | 1413 | 0.89 | -5.97 | 16.99 | -0.12 | 0.33 | 28.72 | 0.78 | 0.43 | 0.72 |
| 3 | summer (JJA) | 1185 | 0.82 | 3.68 | 9.49 | 0.15 | 0.39 | 12.49 | 0.38 | 0.00 | 0.50 |
| 1 | autumn (SON) | 1276 | 0.87 | -1.00 | 15.36 | -0.02 | 0.35 | 21.96 | 0.61 | 0.28 | 0.64 |
| 4 | winter (DJF) | 1945 | 0.81 | 5.96 | 21.21 | 0.11 | 0.39 | 31.17 | 0.72 | 0.35 | 0.68 |
# 12 zmiennych
modStats(mydata = dane3_rf_test,
mod = "mod",
obs = "obs",
type = "season") %>%
knitr::kable(digits = 2)| season | n | FAC2 | MB | MGE | NMB | NMGE | RMSE | r | COE | IOA | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | spring (MAM) | 1413 | 0.89 | -5.64 | 17.15 | -0.11 | 0.33 | 28.76 | 0.77 | 0.43 | 0.71 |
| 3 | summer (JJA) | 1185 | 0.83 | 3.08 | 9.20 | 0.13 | 0.38 | 12.16 | 0.38 | 0.03 | 0.51 |
| 1 | autumn (SON) | 1276 | 0.87 | -0.74 | 15.37 | -0.02 | 0.35 | 21.81 | 0.62 | 0.28 | 0.64 |
| 4 | winter (DJF) | 1945 | 0.80 | 7.07 | 21.68 | 0.13 | 0.40 | 31.64 | 0.71 | 0.34 | 0.67 |
# 13 zmiennych
modStats(mydata = dane4_rf_test,
mod = "mod",
obs = "obs",
type = "season") %>%
knitr::kable(digits = 2)| season | n | FAC2 | MB | MGE | NMB | NMGE | RMSE | r | COE | IOA | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | spring (MAM) | 1413 | 0.89 | -5.82 | 16.94 | -0.11 | 0.33 | 28.68 | 0.78 | 0.44 | 0.72 |
| 3 | summer (JJA) | 1185 | 0.85 | 2.04 | 8.74 | 0.08 | 0.36 | 11.58 | 0.42 | 0.08 | 0.54 |
| 1 | autumn (SON) | 1276 | 0.87 | -0.93 | 15.35 | -0.02 | 0.35 | 21.81 | 0.62 | 0.28 | 0.64 |
| 4 | winter (DJF) | 1945 | 0.78 | 9.23 | 22.86 | 0.17 | 0.42 | 32.24 | 0.71 | 0.30 | 0.65 |
Podsumowujac - wykonując to ćwiczenie poszerzyliśmy nasza wiedzę na temat przeprowadzania analizy danych metorologicznych, oraz o bazie danych GIOŚ. Poznaliśmy wiele możliwości skorzystania z tej wiedzy w celu przeprowadzenia jak najdokladniejszych analiz przy pomocy języka R. W przyszłości wykorystując tą wiedzę oraz zdobyte umiejętności będziemy w stanie osiagnąć zamierzone przez nas cele w stosunkowo krótkim czasie. Uważamy, że najważnieszym aspektem naszej pracy było zdobycie wiedzy z zakresu przeprowadzania prognozy stężenia PM10 przy pomocy RStudio.